#ifndef __CTP23D3D__
#define __CTP23D3D__
#include <stdbool.h>
#include "mesh.h"
#include "enumerations.h"
#include "constants.h"
// ***********************************************************************
// (P2)^3 element, conforming, 3D, 3D value
// ***********************************************************************

/*lhc于2022.1.10号更改至新的编号方式*/

static INT C_T_P2_3D_3D_dof[4] = {3,3,0,0};
static INT C_T_P2_3D_3D_Num_Bas = 30;
static INT C_T_P2_3D_3D_Value_Dim = 3;
static INT C_T_P2_3D_3D_Inter_Dim = 3;  //表示插值值函数的维数
static INT C_T_P2_3D_3D_Polydeg = 2;
static bool C_T_P2_3D_3D_IsOnlyDependOnRefCoord = 1;
static INT C_T_P2_3D_3D_Accuracy = 2;
static MAPTYPE C_T_P2_3D_3D_Maptype = Affine;

static void C_T_P2_3D_3D_InterFun(DOUBLE RefCoord[], INT dim, DOUBLE *values)
{
    INT i;
    for(i = 0; i < 90; i++) values[i] = 0.0;
    for(i=0;i<3;i++)
    {
        values[i*4] = 1.0+2.0*RefCoord[0]*RefCoord[0]+4.0*RefCoord[0]*RefCoord[1]+4.0*RefCoord[0]*RefCoord[2]
        -3*RefCoord[0]+2.0*RefCoord[1]*RefCoord[1]+4.0*RefCoord[1]*RefCoord[2]
        -3*RefCoord[1]+2.0*RefCoord[2]*RefCoord[2]-3.0*RefCoord[2]; 
        values[9+i*4] = -RefCoord[0]+2*RefCoord[0]*RefCoord[0];
        values[18+i*4] = -RefCoord[1]+2*RefCoord[1]*RefCoord[1];
        values[27+i*4] = -RefCoord[2]+2*RefCoord[2]*RefCoord[2];
        values[36+i*4] = 4*RefCoord[0]-4*RefCoord[0]*RefCoord[2]-4*RefCoord[0]*RefCoord[1]-4*RefCoord[0]*RefCoord[0];
        values[45+i*4] = 4*RefCoord[1]-4*RefCoord[1]*RefCoord[2]-4*RefCoord[0]*RefCoord[1]-4*RefCoord[1]*RefCoord[1];
        values[54+i*4] = 4*RefCoord[2]-4*RefCoord[2]*RefCoord[2]-4*RefCoord[0]*RefCoord[2]-4*RefCoord[1]*RefCoord[2];
        values[63+i*4] = 4*RefCoord[0]*RefCoord[1];
        values[72+i*4] = 4*RefCoord[0]*RefCoord[2];
        values[81+i*4] = 4*RefCoord[1]*RefCoord[2];
    }
    // for(i=0;i<3;i++)
    // {
    //     values[i*31] = 1.0+2.0*RefCoord[0]*RefCoord[0]+4.0*RefCoord[0]*RefCoord[1]+4.0*RefCoord[0]*RefCoord[2]
    //     -3*RefCoord[0]+2.0*RefCoord[1]*RefCoord[1]+4.0*RefCoord[1]*RefCoord[2]
    //     -3*RefCoord[1]+2.0*RefCoord[2]*RefCoord[2]-3.0*RefCoord[2]; 
    //     values[i*31+3] = -RefCoord[0]+2*RefCoord[0]*RefCoord[0];
    //     values[i*31+6] = -RefCoord[1]+2*RefCoord[1]*RefCoord[1];
    //     values[i*31+9] = -RefCoord[2]+2*RefCoord[2]*RefCoord[2];
    //     values[i*31+12] = 4*RefCoord[0]-4*RefCoord[0]*RefCoord[2]-4*RefCoord[0]*RefCoord[1]-4*RefCoord[0]*RefCoord[0];
    //     values[i*31+15] = 4*RefCoord[1]-4*RefCoord[1]*RefCoord[2]-4*RefCoord[0]*RefCoord[1]-4*RefCoord[1]*RefCoord[1];
    //     values[i*31+18] = 4*RefCoord[2]-4*RefCoord[2]*RefCoord[2]-4*RefCoord[0]*RefCoord[2]-4*RefCoord[1]*RefCoord[2];
    //     values[i*31+21] = 4*RefCoord[0]*RefCoord[1];
    //     values[i*31+24] = 4*RefCoord[0]*RefCoord[2];
    //     values[i*31+27] = 4*RefCoord[1]*RefCoord[2];
    // }
}


// base function values
static void C_T_P2_3D_3D_D000(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    INT i;
    for(i = 0; i < 90; i++) values[i] = 0.0;
    for(i=0;i<3;i++)
    {
        values[i*4] = 1.0+2.0*RefCoord[0]*RefCoord[0]+4.0*RefCoord[0]*RefCoord[1]+4.0*RefCoord[0]*RefCoord[2]
        -3*RefCoord[0]+2.0*RefCoord[1]*RefCoord[1]+4.0*RefCoord[1]*RefCoord[2]
        -3*RefCoord[1]+2.0*RefCoord[2]*RefCoord[2]-3.0*RefCoord[2]; 
        values[9+i*4] = -RefCoord[0]+2*RefCoord[0]*RefCoord[0];
        values[18+i*4] = -RefCoord[1]+2*RefCoord[1]*RefCoord[1];
        values[27+i*4] = -RefCoord[2]+2*RefCoord[2]*RefCoord[2];
        values[36+i*4] = 4*RefCoord[0]-4*RefCoord[0]*RefCoord[2]-4*RefCoord[0]*RefCoord[1]-4*RefCoord[0]*RefCoord[0];
        values[45+i*4] = 4*RefCoord[1]-4*RefCoord[1]*RefCoord[2]-4*RefCoord[0]*RefCoord[1]-4*RefCoord[1]*RefCoord[1];
        values[54+i*4] = 4*RefCoord[2]-4*RefCoord[2]*RefCoord[2]-4*RefCoord[0]*RefCoord[2]-4*RefCoord[1]*RefCoord[2];
        values[63+i*4] = 4*RefCoord[0]*RefCoord[1];
        values[72+i*4] = 4*RefCoord[0]*RefCoord[2];
        values[81+i*4] = 4*RefCoord[1]*RefCoord[2];
    }
}

// values of the derivatives in RefCoord[0] direction
static void C_T_P2_3D_3D_D100(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    INT i;
    for(i = 0; i < 90; i++) values[i] = 0.0;
    for(i=0;i<3;i++)
    {
        values[i*4] = -3.0+4.0*RefCoord[0]+4.0*RefCoord[1]+4.0*RefCoord[2]; 
        values[9+i*4] = -1.0+4.0*RefCoord[0];
        values[36+i*4] =  4.0-4.0*RefCoord[2]-4.0*RefCoord[1]-8.0*RefCoord[0];
        values[45+i*4] = -4.0*RefCoord[1];
        values[54+i*4] = -4.0*RefCoord[2];
        values[63+i*4] =  4.0*RefCoord[1];
        values[72+i*4] =  4.0*RefCoord[2];
    }
}

// values of the derivatives in RefCoord[1] direction
static void C_T_P2_3D_3D_D010(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    INT i;
    for(i = 0; i < 90; i++) values[i] = 0.0;
    for(i=0;i<3;i++)
    {
        values[i*4] = -3.0+4.0*RefCoord[0]+4.0*RefCoord[1]+4.0*RefCoord[2]; 
        values[18+i*4] = -1.0+4.0*RefCoord[1];
        values[36+i*4] = -4.0*RefCoord[0];
        values[45+i*4] =  4.0-4.0*RefCoord[2]-4.0*RefCoord[0]-8.0*RefCoord[1];
        values[54+i*4] = -4.0*RefCoord[2];
        values[63+i*4] =  4.0*RefCoord[0];
        values[81+i*4] =  4.0*RefCoord[2];
    }
}
// values of the derivatives in RefCoord[2] direction
static void C_T_P2_3D_3D_D001(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    INT i;
    for(i = 0; i < 90; i++) values[i] = 0.0;
    for(i=0;i<3;i++)
    {
        values[i*4] = -3.0+4.0*RefCoord[0]+4.0*RefCoord[1]+4.0*RefCoord[2]; 
        values[27+i*4] = -1.0+4.0*RefCoord[2];
        values[36+i*4] = -4.0*RefCoord[0];
        values[45+i*4] = -4.0*RefCoord[1];
        values[54+i*4] =  4.0-4.0*RefCoord[0]-4.0*RefCoord[1]-8.0*RefCoord[2];
        values[72+i*4] =  4.0*RefCoord[0];
        values[81+i*4] =  4.0*RefCoord[1];
    }
}

// values of the derivatives in RefCoord[0]-RefCoord[0]  direction
static void C_T_P2_3D_3D_D200(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    INT i; 
    for(i = 0; i < 90; i++) values[i] = 0.0;

    for(i=0;i<3;i++)
    {
        values[i*4] =  4.0; 
        values[9+i*4] =  4.0;
        values[36+i*4] = -8.0;
}
}
// values of the derivatives in RefCoord[0]-RefCoord[0]  direction
static void C_T_P2_3D_3D_D020(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    INT i; 
    for(i = 0; i < 90; i++) values[i] = 0.0;
    for(i=0;i<3;i++)
    {  
        values[i*4] =  4.0; 
        values[18+i*4] =  4.0;
        values[45+i*4] = -8.0;
    }
}
// values of the derivatives in RefCoord[0]-RefCoord[0]  direction
static void C_T_P2_3D_3D_D002(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    INT i;
    for(i = 0; i < 90; i++) values[i] = 0.0;
    for(i=0;i<3;i++)
    {
        values[i*4] =  4.0; 
        values[27+i*4] =  4.0;
        values[54+i*4] = -8.0;
    }
}
// values of the derivatives in RefCoord[0]-RefCoord[0]  direction
static void C_T_P2_3D_3D_D110(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    INT i; 
    for(i = 0; i < 90; i++) values[i] = 0.0;
    for(i=0;i<3;i++)
    {
        values[i*4] =  4.0; 
        values[36+i*4] = -4.0;
        values[45+i*4] = -4.0;
        values[63+i*4] =  4.0;
    }
}
// values of the derivatives in RefCoord[0]-RefCoord[0]  direction
static void C_T_P2_3D_3D_D101(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    INT i; 
    for(i = 0; i < 90; i++) values[i] = 0.0;
    for(i=0;i<3;i++)
    {
        values[i*4] =  4.0; 
        values[36+i*4] = -4.0;
        values[54+i*4] = -4.0;
        values[72+i*4] =  4.0;
    }
}
// values of the derivatives in RefCoord[0]-RefCoord[0]  direction
static void C_T_P2_3D_3D_D011(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
    INT i; 
    for(i = 0; i < 90; i++) values[i] = 0.0;
    for(i=0;i<3;i++)
    {
        values[i*4] =  4.0; 
        values[45+i*4] = -4.0;
        values[54+i*4] = -4.0;
        values[81+i*4] =  4.0;
    }
}
static void C_T_P2_3D_3D_Nodal(ELEMENT * elem, FUNCTIONVEC *fun, INT dim, DOUBLE* values)
{
    //注意fun是一个向量函数
    DOUBLE coord[3];
    INT i;
    //先计算4个节点上的函数值
    for(i=0;i<4;i++)
    {
       // dof on the verts 
        coord[0] = elem->Vert_X[i]; 
        coord[1] = elem->Vert_Y[i];
        coord[2] = elem->Vert_Z[i];
        fun(coord, dim, values+i*dim);
    }
    INT vert0[6] = {0, 0, 0, 1, 1, 2};
    INT vert1[6] = {1, 2, 3, 2, 3, 3};

    for(i=0;i<6;i++)
    {
       coord[0] = 0.5*(elem->Vert_X[vert0[i]]+elem->Vert_X[vert1[i]]); 
       coord[1] = 0.5*(elem->Vert_Y[vert0[i]]+elem->Vert_Y[vert1[i]]); 
       coord[2] = 0.5*(elem->Vert_Z[vert0[i]]+elem->Vert_Z[vert1[i]]); 
       fun(coord, dim, values+(4+i)*dim);
    } 
}
#endif
